//
// Created by baihe on 2016-07-01
// Update on 2019-07
//

#include "slave.h"
#include <math.h>
#include <time.h>
  
__gsm__ volatile unsigned long reply, _reply[2], put_reply;
__gsm__ int my_id;
__gsm__ int Nvalue;
__gsm__ int bigblocknum;
__gsm__ int xstart, ystart, zstart;
__gsm__ int nlocalx, nlocaly, nlocalz;
__gsm__ int nghostx, nghosty;
__gsm__ int NeighborOffset[68];
__gsm__ double _spline_[5005];
__gsm__ int block_size;
__gsm__ int nx_per_step = 16;
__gsm__ int ny_per_step = 5;
__gsm__ int nz_per_step = 3;
__gsm__ int natoms_per_step = 240;  //nx * ny * nz
__gsm__ int nxy_per_step = 80;  //nx * ny
__gsm__ int nx_per_step_cal = 8;

#define max(x, y) (x>y?x:y)
#define min(x, y) (x>y?y:x)

struct AtomElement {
    unsigned long id;
    enum {
        INVALID = -1, Fe = 0
    } type;
    double x[3];
    double v[3];
    double f[3];
    double rho;
    double df;
};

void rtc(unsigned long *counter);

int BeginOfArea(int natom, int start, int subindex);

int EndOfArea(int natom, int start, int subindex);

long int IndexOf3DIndex(long int xIndex, long int yIndex, long int zIndex);

long int co_IndexOf3DIndex(long int xIndex, long int yIndex, long int zIndex);

void calculateNeighbourIndices();

void init(int *a) {
    int i, j;
    my_id = get_thread_id();
    xstart = a[0];
    ystart = a[1];
    zstart = a[2];
    nlocalx = a[3];
    nlocaly = a[4];
    nlocalz = a[5];
    nghostx = a[6];
    nghosty = a[7];
    Nvalue = a[8] * 2;
    block_size = nlocalz / Nvalue;
    bigblocknum = (nlocalz - Nvalue * block_size) ? 1 : 0;
    calculateNeighbourIndices();

}

void calculateNeighbourIndices() {
    int i;
    double x, y, z;
    int xIndex, yIndex, zIndex;
    int index = 0;
    int mark = 0;
    double _cutoffRadius = 5.6;
    double _latticeconst = 2.85532;
    int _cutlattice = 2;
    double cut_times_latti = _cutoffRadius / _latticeconst;
    double square, r;
    long int offset;
    for (zIndex = -_cutlattice; zIndex <= _cutlattice; zIndex++) {
        for (yIndex = -_cutlattice; yIndex <= _cutlattice; yIndex++) {
            for (xIndex = -_cutlattice * 2; xIndex <= _cutlattice * 2; xIndex++) {
                z = (double) zIndex + (((double) (xIndex % 2)) / 2);
                y = (double) yIndex + (((double) (xIndex % 2)) / 2);
                x = (double) xIndex / 2;
                square = x * x + y * y + z * z;
                offset;
                r = sqrt(square);
                if (r < (cut_times_latti + 0.4)) {
                    offset = co_IndexOf3DIndex(xIndex, yIndex, zIndex);
                    if (offset > 0) {
                        NeighborOffset[index] = offset;
                        index++;
                    }
                }
                z = (double) zIndex - (((double) (xIndex % 2)) / 2);
                y = (double) yIndex - (((double) (xIndex % 2)) / 2);
                x = (double) xIndex / 2;
                square = x * x + y * y + z * z;
                r = sqrt(square);
                if (r < (cut_times_latti + 0.4)) {
                    offset = co_IndexOf3DIndex(xIndex, yIndex, zIndex);
                    if (offset > 0) {
                        for (i = 0; i < index; i++) {
                            if (NeighborOffset[i] == offset) {
                                mark = 1;
                                break;
                            }
                        }
                        if (mark != 1) {
                            NeighborOffset[index] = offset;
                            index++;
                        }
                        mark = 0;
                    }
                }
            }
        }
    }
}

long int co_IndexOf3DIndex(long int xIndex, long int yIndex, long int zIndex) {
    return (zIndex * ny_per_step + yIndex) * nx_per_step + xIndex;
}

void init_spline(double **spline[]) {
    //rho_spline = spline[0];
    //f_spline = spline[1];
    //phi_spline = spline[2];
}

void cal_rho1(double cutoffRadius,int rho_n,double invDx,
            struct AtomElement *atomlist, double *value) 
{
    int subindex_z, zOmp_start, zOmp_end;
    int i, j, k, l, m_block, kk;
    // double cutoffRadius;
    double xtemp, ytemp, ztemp;
    double delx, dely, delz;
    int n;
    double *list, (*spline)[7], *spline_mem;
    // struct AtomElement *atomlist;
    struct AtomElement atoms[natoms_per_step];
    // double *value;
    double _spline[4];
    // int rho_n;
    // double invDx;
    double dist2, r, p, rhoTmp;
    int m;
    int m_step, x_end;
    int n_block;
    unsigned long start, stop, time, ts, te;
    unsigned long comms, comme, commt;
    unsigned int ch,ch1,ch2,ch3,ch4,ch5,ch6;

    int ntime = 0;
    commt = 0;
    time = 0;

    // list = a[0];
    // cutoffRadius = *a[0];
    // rho_n = *a[2];
    // invDx = *a[3];
    // value = a[4];

    // atomlist = (struct AtomElement *) list;
    subindex_z = (my_id * 2) % Nvalue;
    zOmp_start = BeginOfArea(nlocalz, zstart, subindex_z);
    zOmp_end = EndOfArea(nlocalz, zstart, subindex_z);
    m_step = 0;
    n_block = (nlocalx % nx_per_step_cal) ? (nlocalx / nx_per_step_cal + 1) : (nlocalx / nx_per_step_cal);
    rtc(&ts);
    // reply = 0;
    // athread_get(PE_MODE, &value[0], &_spline_[0], 5000 * 8, (void *) &reply, 0, 0, 0);
    // while (reply != 1);
	ch = dma_p2p(&value[0], 1, 5000*sizeof(long), 0,
		 &_spline_[0], 1, 5000*sizeof(long), 0,
		false, 0);
    dma_wait(ch);
    for (k = zOmp_start; k < zOmp_end; k++) {
        for (j = ystart; j < nlocaly + ystart; j++) {
            for (m_step = 0; m_step < n_block; m_step++) {
                rtc(&comms);
                // reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[0], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[0], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step * 2], nxy_per_step * 104, (void *) &reply,
                //             0, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = nx_per_step * 2 + xstart - 2;
                if ((m_step % n_block) < n_block - 1)
                    x_end = nx_per_step_cal;
                else
                    x_end = nlocalx - nx_per_step_cal * (n_block - 1);
                // while (reply != 3);
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
                rtc(&comme);
                commt += comme - comms;
                for (i = 0; i < x_end; i++) {
                    if (atoms[kk].type != -1) {
                        xtemp = atoms[kk].x[0];
                        ytemp = atoms[kk].x[1];
                        ztemp = atoms[kk].x[2];
                        for (l = 0; l < 68; l++) {
                            n = kk + NeighborOffset[l];
                            if (atoms[n].type == -1) continue;
                            delx = xtemp - atoms[n].x[0];
                            dely = ytemp - atoms[n].x[1];
                            delz = ztemp - atoms[n].x[2];
                            dist2 = delx * delx + dely * dely + delz * delz;
                            if (dist2 < (cutoffRadius * cutoffRadius)) {
                                rtc(&start);
                                r = sqrt(dist2);
                                p = r * invDx + 1.0;
                                m = (int) p;
                                m = max(1, min(m, (rho_n - 1)));

                                p -= m;
                                p = min(p, 1.0);
                                _spline[3] = _spline_[m];
                                _spline[2] = ((_spline_[m - 2] - _spline_[m + 2]) +
                                              8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0;
                                _spline[1] = 3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * _spline[2] -
                                             ((_spline_[m - 1] - _spline_[m + 3]) +
                                              8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0;
                                _spline[0] = _spline[2] + ((_spline_[m - 1] - _spline_[m + 3]) +
                                                           8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0 -
                                             2.0 * (_spline_[m + 1] - _spline_[m]);
                                rtc(&stop);
                                time += stop - start;
                                ntime++;
                                rhoTmp = ((_spline[0] * p + _spline[1]) * p + _spline[2]) * p + _spline[3];

                                atoms[kk].rho += rhoTmp;
                                atoms[n].rho += rhoTmp;
                            }
                        }
                    }
                    kk++;
                }
                rtc(&comms);
                put_reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_put(PE_MODE, &atoms[0], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atoms[0], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);    
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_put(PE_MODE, &atoms[nxy_per_step], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);   
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                ch3 = dma_p2p(&atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0); 
                // athread_put(PE_MODE, &atoms[nxy_per_step * 2], &atomlist[kk] + 2, nxy_per_step * 104,
                //             (void *) &put_reply, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                // while (put_reply != 3);
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
                rtc(&comme);
                commt += comme - comms;
            }
        }
    }
    rtc(&te);
    if (my_id == 0) printf("total: %ld, compute: %ld %d, commication: %d\n", te - ts, time, ntime, commt);
}

void cal_rho2(double *a[]) {
    int subindex_z, zOmp_start, zOmp_end;
    int i, j, k, l, kk, m_block;
    double cutoffRadius;
    double xtemp, ytemp, ztemp;
    double delx, dely, delz;
    int n;
    double *list, (*spline)[7], *spline_mem;
    struct AtomElement *atomlist;
    struct AtomElement atoms[natoms_per_step];
    double _spline[4];
    int rho_n;
    double invDx;
    double dist2, r, p, rhoTmp;
    int m;
    int m_step, x_end;
    int n_block;
    unsigned int ch,ch1,ch2,ch3,ch4,ch5,ch6;
    list = a[0];
    cutoffRadius = *a[1];
    rho_n = *a[2];
    invDx = *a[3];

    atomlist = (struct AtomElement *) list;
    subindex_z = (my_id * 2 + 1) % Nvalue;
    zOmp_start = BeginOfArea(nlocalz, zstart, subindex_z);
    zOmp_end = EndOfArea(nlocalz, zstart, subindex_z);
    m_step = 0;
    n_block = (nlocalx % nx_per_step_cal) ? (nlocalx / nx_per_step_cal + 1) : (nlocalx / nx_per_step_cal);


    for (k = zOmp_start; k < zOmp_end; k++) {
        for (j = ystart; j < nlocaly + ystart; j++) {
            for (m_step = 0; m_step < n_block; m_step++) {
                reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[0], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[0], 1, nxy_per_step * 104, 0,
                             false, 0);
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step * 2], nxy_per_step * 104, (void *) &reply,
                //             0, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = nx_per_step * 2 + xstart - 2;
                if ((m_step % n_block) < n_block - 1)
                    x_end = nx_per_step_cal;
                else
                    x_end = nlocalx - nx_per_step_cal * (n_block - 1);
                // while (reply != 3);
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
                for (i = 0; i < x_end; i++) {
                    if (atoms[kk].type != -1) {
                        xtemp = atoms[kk].x[0];
                        ytemp = atoms[kk].x[1];
                        ztemp = atoms[kk].x[2];
                        for (l = 0; l < 68; l++) {
                            n = kk + NeighborOffset[l];
                            if (atoms[n].type == -1) continue;
                            delx = xtemp - atoms[n].x[0];
                            dely = ytemp - atoms[n].x[1];
                            delz = ztemp - atoms[n].x[2];
                            dist2 = delx * delx + dely * dely + delz * delz;
                            if (dist2 < (cutoffRadius * cutoffRadius)) {
                                r = sqrt(dist2);
                                p = r * invDx + 1.0;
                                m = (int) p;
                                m = max(1, min(m, (rho_n - 1)));
                                p -= m;
                                p = min(p, 1.0);
                                _spline[3] = _spline_[m];
                                _spline[2] = ((_spline_[m - 2] - _spline_[m + 2]) +
                                              8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0;
                                _spline[1] = 3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * _spline[2] -
                                             ((_spline_[m - 1] - _spline_[m + 3]) +
                                              8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0;
                                _spline[0] = _spline[2] + ((_spline_[m - 1] - _spline_[m + 3]) +
                                                           8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0 -
                                             2.0 * (_spline_[m + 1] - _spline_[m]);
                                rhoTmp = ((_spline[0] * p + _spline[1]) * p + _spline[2]) * p + _spline[3];
                                atoms[kk].rho += rhoTmp;
                                atoms[n].rho += rhoTmp;
                            }
                        }
                    }
                    kk++;
                }
                put_reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_put(PE_MODE, &atoms[0], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atoms[0], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);                   
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_put(PE_MODE, &atoms[nxy_per_step], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);                   
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);   
                // athread_put(PE_MODE, &atoms[nxy_per_step * 2], &atomlist[kk] + 2, nxy_per_step * 104,
                //             (void *) &put_reply, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                // while (put_reply != 3);
                ch3 = dma_p2p(&atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0); 
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);            
            }
        }
    }
}

void cal_df(double *a[]) {
    int zOmp_start, zOmp_end;
    int i, j, k, kk;
    double cutoffRadius;
    double *list, (*spline)[7];
    double *value;
    int df_n;
    double invDx;
    int m;
    double dist2, p, dfEmbed;
    struct AtomElement *atomlist;
    struct AtomElement atoms[natoms_per_step];
    double _spline[3];
    int m_step, x_end;
    int n_block;
    int natoms;
    unsigned int ch,ch1,ch2,ch3,ch4,ch5,ch6;
    list = a[0];
    cutoffRadius = *a[1];
    df_n = *a[2];
    invDx = *a[3];
    value = a[4];

    atomlist = (struct AtomElement *) list;
    zOmp_start = zstart + ((nlocalz / (Nvalue / 2)) * my_id);
    zOmp_end = zstart + ((nlocalz / (Nvalue / 2)) * (my_id + 1));
    if (my_id == (Nvalue / 2 - 1)) {
        zOmp_end = nlocalz + zstart;
    }
    n_block = (nlocalx % natoms_per_step) ? (nlocalx / natoms_per_step + 1) : (nlocalx / natoms_per_step);

    // reply = 0;
    // athread_get(PE_MODE, &value[0], &_spline_[0], 5000 * 8, (void *) &reply, 0, 0, 0);
    // while (reply != 1);
	ch = dma_p2p(&value[0], 1, 5000*sizeof(long), 0,
		 &_spline_[0], 1, 5000*sizeof(long), 0,
		false, 0);
    dma_wait(ch);
    for (k = zOmp_start; k < zOmp_end; k++) {
        for (j = ystart; j < nlocaly + ystart; j++) {
            for (m_step = 0; m_step < n_block; m_step++) {
                reply = 0;
                kk = IndexOf3DIndex(xstart + m_step * natoms_per_step, j, k);
                if (nlocalx >= (m_step + 1) * natoms_per_step) natoms = natoms_per_step;
                else natoms = nlocalx - (m_step) * natoms_per_step;
                //athread_get(PE_MODE, &atomlist[kk], &atoms[0], natoms * 104, (void *) &reply, 0, 0, 0);
                ch1 = dma_p2p(&atomlist[kk], 1,  natoms * 104, 0 ,
                             &atoms[0], 1,  natoms * 104, 0,
                             false, 0); 
                // while (reply != 1);
                dma_wait(ch1);
                for (i = 0; i < natoms; i++) {
                    if (atoms[i].type == -1) continue;
                    p = atoms[i].rho * invDx + 1.0;
                    m = (int) p;
                    m = max(1, min(m, (df_n - 1)));
                    p -= m;
                    p = min(p, 1.0);
                    _spline[2] =
                            (((_spline_[m - 2] - _spline_[m + 2]) + 8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) *
                            invDx;
                    _spline[1] =
                            2.0 * (3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * (((_spline_[m - 2] - _spline_[m + 2]) +
                                                                                   8.0 * (_spline_[m + 1] -
                                                                                          _spline_[m - 1])) / 12.0) -
                                   (((_spline_[m - 1] - _spline_[m + 3]) +
                                     8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0)) * invDx;
                    _spline[0] =
                            3.0 * ((((_spline_[m - 2] - _spline_[m + 2]) + 8.0 * (_spline_[m + 1] - _spline_[m - 1])) /
                                    12.0) +
                                   (((_spline_[m - 1] - _spline_[m + 3]) + 8.0 * (_spline_[m + 2] - _spline_[m])) /
                                    12.0) - 2.0 * (_spline_[m + 1] - _spline_[m])) * invDx;

                    dfEmbed = (_spline[0] * p + _spline[1]) * p + _spline[2];

                    atoms[i].df = dfEmbed;
                }
                put_reply = 0;
                kk = IndexOf3DIndex(xstart + m_step * natoms_per_step, j, k);
                //athread_put(PE_MODE, &atoms[0], &atomlist[kk], natoms * 104, (void *) &put_reply, 0, 0);
                // while (put_reply != 1);
                ch1 = dma_p2p(&atoms[0], 1,  natoms * 104, 0 ,
                             &atomlist[kk], 1,  natoms * 104, 0,
                             false, 0); 
                dma_wait(ch1);
            }
        }
    }
}

void cal_force1(double *a[]) {
    int subindex_z, zOmp_start, zOmp_end;
    int i, j, k, l, m_block, kk;
    double cutoffRadius;
    double xtemp, ytemp, ztemp;
    double delx, dely, delz;
    int n;
    double *list, (*spline)[7];
    struct AtomElement *atomlist;
    struct AtomElement atoms[natoms_per_step];
    double *value;
    double _spline[7];
    int phi_n;
    double invDx;
    double dist2, r, p, phiTmp, dPhi, dRho;
    double z2, z2p, recip, phi, phip, psip, fpair;
    int m;
    int m_step, x_end;
    int n_block;
    unsigned int ch,ch1,ch2,ch3,ch4,ch5,ch6;
    list = a[0];
    cutoffRadius = *a[1];
    phi_n = *a[2];
    invDx = *a[3];
    value = a[4];

    atomlist = (struct AtomElement *) list;
    subindex_z = (my_id * 2) % Nvalue;
    zOmp_start = BeginOfArea(nlocalz, zstart, subindex_z);
    zOmp_end = EndOfArea(nlocalz, zstart, subindex_z);
    m_step = 0;
    n_block = (nlocalx % nx_per_step_cal) ? (nlocalx / nx_per_step_cal + 1) : (nlocalx / nx_per_step_cal);

    reply = 0;
    // athread_get(PE_MODE, &value[0], &_spline_[0], 5000 * 8, (void *) &reply, 0, 0, 0);
    // while (reply != 1);
	ch = dma_p2p(&value[0], 1, 5000*sizeof(long), 0,
		 &_spline_[0], 1, 5000*sizeof(long), 0,
		false, 0);
    dma_wait(ch);
    for (k = zOmp_start; k < zOmp_end; k++) {
        for (j = ystart; j < nlocaly + ystart; j++) {
            for (m_step = 0; m_step < n_block; m_step++) {
                reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[0], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[0], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step * 2], nxy_per_step * 104, (void *) &reply,
                //             0, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = nx_per_step * 2 + xstart - 2;
                if ((m_step % n_block) < n_block - 1)
                    x_end = nx_per_step_cal;
                else
                    x_end = nlocalx - nx_per_step_cal * (n_block - 1);
                // while (reply != 3);
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
                for (i = 0; i < x_end; i++) {
                    if (atoms[kk].type != -1) {
                        xtemp = atoms[kk].x[0];
                        ytemp = atoms[kk].x[1];
                        ztemp = atoms[kk].x[2];
                        for (l = 0; l < 68; l++) {
                            n = kk + NeighborOffset[l];
                            if (atoms[n].type == -1) continue;
                            delx = xtemp - atoms[n].x[0];
                            dely = ytemp - atoms[n].x[1];
                            delz = ztemp - atoms[n].x[2];
                            dist2 = delx * delx + dely * dely + delz * delz;
                            //printf("nei:%lf, %lf, %lf\n", _x[n * 3], _x[n * 3+2], _x[n * 3+2]);
                            if (dist2 < (cutoffRadius * cutoffRadius)) {
                                r = sqrt(dist2);
                                p = r * invDx + 1.0;
                                m = (int) p;
                                m = max(1, min(m, (phi_n - 1)));

                                p -= m;
                                p = min(p, 1.0);
                                //while(reply!=1);
                                _spline[6] = _spline_[m];
                                _spline[5] = ((_spline_[m - 2] - _spline_[m + 2]) +
                                              8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0;
                                _spline[4] = 3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * _spline[5] -
                                             ((_spline_[m - 1] - _spline_[m + 3]) +
                                              8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0;
                                _spline[3] = _spline[5] + ((_spline_[m - 1] - _spline_[m + 3]) +
                                                           8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0 -
                                             2.0 * (_spline_[m + 1] - _spline_[m]);
                                _spline[2] = (((_spline_[m - 2] - _spline_[m + 2]) +
                                               8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) * invDx;
                                _spline[1] = 2.0 * (3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * (((_spline_[m - 2] -
                                                                                                     _spline_[m + 2]) +
                                                                                                    8.0 *
                                                                                                    (_spline_[m + 1] -
                                                                                                     _spline_[m - 1])) /
                                                                                                   12.0) -
                                                    (((_spline_[m - 1] - _spline_[m + 3]) +
                                                      8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0)) * invDx;
                                _spline[0] = 3.0 * ((((_spline_[m - 2] - _spline_[m + 2]) +
                                                      8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) +
                                                    (((_spline_[m - 1] - _spline_[m + 3]) +
                                                      8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0) -
                                                    2.0 * (_spline_[m + 1] - _spline_[m])) * invDx;
                                phiTmp = ((_spline[3] * p + _spline[4]) * p + _spline[5]) * p + _spline[6];
                                dPhi = (_spline[0] * p + _spline[1]) * p + _spline[2];

                                z2 = phiTmp;
                                z2p = dPhi;
                                recip = 1.0 / r;
                                phi = z2 * recip;
                                phip = z2p * recip - phi * recip;
                                psip = phip;
                                fpair = -psip * recip;

                                atoms[kk].f[0] += delx * fpair;
                                atoms[kk].f[1] += dely * fpair;
                                atoms[kk].f[2] += delz * fpair;
                                atoms[n].f[0] -= delx * fpair;
                                atoms[n].f[1] -= dely * fpair;
                                atoms[n].f[2] -= delz * fpair;
                            }
                        }
                    }
                    kk++;
                }
                put_reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_put(PE_MODE, &atoms[0], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atoms[0], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);    
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_put(PE_MODE, &atoms[nxy_per_step], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);   
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_put(PE_MODE, &atoms[nxy_per_step * 2], &atomlist[kk] + 2, nxy_per_step * 104,
                //             (void *) &put_reply, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0); 
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
                // while (put_reply != 3);
            }
        }
    }
}

void cal_force2(double *a[]) {
    int subindex_z, zOmp_start, zOmp_end;
    int i, j, k, l, kk, m_block;
    double cutoffRadius;
    double xtemp, ytemp, ztemp;
    double delx, dely, delz;
    int n;
    double *list, (*spline)[7];
    struct AtomElement *atomlist;
    struct AtomElement atoms[natoms_per_step];
    double *value;
    double _spline[7];
    int phi_n;
    double invDx;
    double dist2, r, p, phiTmp, dPhi, dRho;
    double z2, z2p, recip, phi, phip, psip, fpair;
    int m;
    int m_step, x_end;
    int n_block;
    unsigned int ch,ch1,ch2,ch3,ch4,ch5,ch6;
    list = a[0];
    cutoffRadius = *a[1];
    phi_n = *a[2];
    invDx = *a[3];
    value = a[4];

    atomlist = (struct AtomElement *) list;
    subindex_z = (my_id * 2 + 1) % Nvalue;
    zOmp_start = BeginOfArea(nlocalz, zstart, subindex_z);
    zOmp_end = EndOfArea(nlocalz, zstart, subindex_z);
    m_step = 0;
    n_block = (nlocalx % nx_per_step_cal) ? (nlocalx / nx_per_step_cal + 1) : (nlocalx / nx_per_step_cal);

    for (k = zOmp_start; k < zOmp_end; k++) {
        for (j = ystart; j < nlocaly + ystart; j++) {
            for (m_step = 0; m_step < n_block; m_step++) {
                reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[0], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[0], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step * 2], nxy_per_step * 104, (void *) &reply,
                //             0, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = nx_per_step * 2 + xstart - 2;
                if ((m_step % n_block) < n_block - 1)
                    x_end = nx_per_step_cal;
                else
                    x_end = nlocalx - nx_per_step_cal * (n_block - 1);
                // while (reply != 3);
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
                for (i = 0; i < x_end; i++) {
                    if (atoms[kk].type != -1) {
                        xtemp = atoms[kk].x[0];
                        ytemp = atoms[kk].x[1];
                        ztemp = atoms[kk].x[2];
                        for (l = 0; l < 68; l++) {
                            n = kk + NeighborOffset[l];
                            if (atoms[n].type == -1) continue;
                            delx = xtemp - atoms[n].x[0];
                            dely = ytemp - atoms[n].x[1];
                            delz = ztemp - atoms[n].x[2];
                            dist2 = delx * delx + dely * dely + delz * delz;
                            if (dist2 < (cutoffRadius * cutoffRadius)) {
                                r = sqrt(dist2);
                                p = r * invDx + 1.0;
                                m = (int) p;
                                m = max(1, min(m, (phi_n - 1)));

                                p -= m;
                                p = min(p, 1.0);
                                _spline[6] = _spline_[m];
                                _spline[5] = ((_spline_[m - 2] - _spline_[m + 2]) +
                                              8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0;
                                _spline[4] = 3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * _spline[5] -
                                             ((_spline_[m - 1] - _spline_[m + 3]) +
                                              8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0;
                                _spline[3] = _spline[5] + ((_spline_[m - 1] - _spline_[m + 3]) +
                                                           8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0 -
                                             2.0 * (_spline_[m + 1] - _spline_[m]);
                                _spline[2] = (((_spline_[m - 2] - _spline_[m + 2]) +
                                               8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) * invDx;
                                _spline[1] = 2.0 * (3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * (((_spline_[m - 2] -
                                                                                                     _spline_[m + 2]) +
                                                                                                    8.0 *
                                                                                                    (_spline_[m + 1] -
                                                                                                     _spline_[m - 1])) /
                                                                                                   12.0) -
                                                    (((_spline_[m - 1] - _spline_[m + 3]) +
                                                      8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0)) * invDx;
                                _spline[0] = 3.0 * ((((_spline_[m - 2] - _spline_[m + 2]) +
                                                      8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) +
                                                    (((_spline_[m - 1] - _spline_[m + 3]) +
                                                      8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0) -
                                                    2.0 * (_spline_[m + 1] - _spline_[m])) * invDx;
                                phiTmp = ((_spline[3] * p + _spline[4]) * p + _spline[5]) * p + _spline[6];
                                dPhi = (_spline[0] * p + _spline[1]) * p + _spline[2];

                                z2 = phiTmp;
                                z2p = dPhi;
                                recip = 1.0 / r;
                                phi = z2 * recip;
                                phip = z2p * recip - phi * recip;
                                psip = phip;
                                fpair = -psip * recip;

                                atoms[kk].f[0] += delx * fpair;
                                atoms[kk].f[1] += dely * fpair;
                                atoms[kk].f[2] += delz * fpair;
                                atoms[n].f[0] -= delx * fpair;
                                atoms[n].f[1] -= dely * fpair;
                                atoms[n].f[2] -= delz * fpair;
                            }
                        }
                    }
                    kk++;
                }
                put_reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_put(PE_MODE, &atoms[0], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atoms[0], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);    
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_put(PE_MODE, &atoms[nxy_per_step], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);  
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_put(PE_MODE, &atoms[nxy_per_step * 2], &atomlist[kk] + 2, nxy_per_step * 104,
                //             (void *) &put_reply, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0); 
                // while (put_reply != 3);
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
            }
        }
    }

    value = a[5];
    // reply = 0;
    // athread_get(PE_MODE, &value[0], &_spline_[0], 5000 * 8, (void *) &reply, 0, 0, 0);
    // while (reply != 1);
	ch = dma_p2p(&value[0], 1, 5000*sizeof(long), 0,
		 &_spline_[0], 1, 5000*sizeof(long), 0,
		false, 0);
    dma_wait(ch);
    for (k = zOmp_start; k < zOmp_end; k++) {
        for (j = ystart; j < nlocaly + ystart; j++) {
            for (m_step = 0; m_step < n_block; m_step++) {
                reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[0], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[0], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step * 2], nxy_per_step * 104, (void *) &reply,
                //             0, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = nx_per_step * 2 + xstart - 2;
                if ((m_step % n_block) < n_block - 1)
                    x_end = nx_per_step_cal;
                else
                    x_end = nlocalx - nx_per_step_cal * (n_block - 1);
                // while (reply != 3);
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
                for (i = 0; i < x_end; i++) {
                    if (atoms[kk].type != -1) {
                        xtemp = atoms[kk].x[0];
                        ytemp = atoms[kk].x[1];
                        ztemp = atoms[kk].x[2];
                        for (l = 0; l < 68; l++) {
                            n = kk + NeighborOffset[l];
                            if (atoms[n].type == -1) continue;
                            delx = xtemp - atoms[n].x[0];
                            dely = ytemp - atoms[n].x[1];
                            delz = ztemp - atoms[n].x[2];
                            dist2 = delx * delx + dely * dely + delz * delz;
                            if (dist2 < (cutoffRadius * cutoffRadius)) {
                                r = sqrt(dist2);
                                p = r * invDx + 1.0;
                                m = (int) p;
                                m = max(1, min(m, (phi_n - 1)));
                                p -= m;
                                p = min(p, 1.0);
                                _spline[2] = (((_spline_[m - 2] - _spline_[m + 2]) +
                                               8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) * invDx;
                                _spline[1] = 2.0 * (3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * (((_spline_[m - 2] -
                                                                                                     _spline_[m + 2]) +
                                                                                                    8.0 *
                                                                                                    (_spline_[m + 1] -
                                                                                                     _spline_[m - 1])) /
                                                                                                   12.0) -
                                                    (((_spline_[m - 1] - _spline_[m + 3]) +
                                                      8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0)) * invDx;
                                _spline[0] = 3.0 * ((((_spline_[m - 2] - _spline_[m + 2]) +
                                                      8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) +
                                                    (((_spline_[m - 1] - _spline_[m + 3]) +
                                                      8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0) -
                                                    2.0 * (_spline_[m + 1] - _spline_[m])) * invDx;

                                dRho = (_spline[0] * p + _spline[1]) * p + _spline[2];
                                recip = 1.0 / r;
                                psip = (atoms[kk].df + atoms[n].df) * dRho;
                                fpair = -psip * recip;

                                atoms[kk].f[0] += delx * fpair;
                                atoms[kk].f[1] += dely * fpair;
                                atoms[kk].f[2] += delz * fpair;
                                atoms[n].f[0] -= delx * fpair;
                                atoms[n].f[1] -= dely * fpair;
                                atoms[n].f[2] -= delz * fpair;
                            }
                        }
                    }
                    kk++;
                }
                put_reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_put(PE_MODE, &atoms[0], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atoms[0], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);                   
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_put(PE_MODE, &atoms[nxy_per_step], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);                  
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_put(PE_MODE, &atoms[nxy_per_step * 2], &atomlist[kk] + 2, nxy_per_step * 104,
                //             (void *) &put_reply, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);         
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);       
                // while (put_reply != 3);
            }
        }
    }
}

void cal_force3(double *a[]) {
    int subindex_z, zOmp_start, zOmp_end;
    int i, j, k, l, m_block, kk;
    double cutoffRadius;
    double xtemp, ytemp, ztemp;
    double delx, dely, delz;
    int n;
    double *list, (*spline)[7];
    struct AtomElement *atomlist;
    struct AtomElement atoms[natoms_per_step];
    double *value;
    double _spline[7];
    int phi_n;
    double invDx;
    double dist2, r, p, phiTmp, dPhi, dRho;
    double z2, z2p, recip, phi, phip, psip, fpair;
    int m;
    int m_step, x_end;
    int n_block;
    unsigned int ch,ch1,ch2,ch3,ch4,ch5,ch6;
    list = a[0];
    cutoffRadius = *a[1];
    phi_n = *a[2];
    invDx = *a[3];

    atomlist = (struct AtomElement *) list;
    subindex_z = (my_id * 2) % Nvalue;
    zOmp_start = BeginOfArea(nlocalz, zstart, subindex_z);
    zOmp_end = EndOfArea(nlocalz, zstart, subindex_z);
    m_step = 0;
    n_block = (nlocalx % nx_per_step_cal) ? (nlocalx / nx_per_step_cal + 1) : (nlocalx / nx_per_step_cal);

    for (k = zOmp_start; k < zOmp_end; k++) {
        for (j = ystart; j < nlocaly + ystart; j++) {
            for (m_step = 0; m_step < n_block; m_step++) {
                reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[0], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[0], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step], nxy_per_step * 104, (void *) &reply, 0,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                athread_get(PE_MODE, &atomlist[kk] + 2, &atoms[nxy_per_step * 2], nxy_per_step * 104, (void *) &reply,
                            0, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atomlist[kk] + 2, ny_per_step,  nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             &atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                             false, 0); 
                kk = nx_per_step * 2 + xstart - 2;
                if ((m_step % n_block) < n_block - 1)
                    x_end = nx_per_step_cal;
                else
                    x_end = nlocalx - nx_per_step_cal * (n_block - 1);
                // while (reply != 3);
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);
                for (i = 0; i < x_end; i++) {
                    if (atoms[kk].type != -1) {
                        xtemp = atoms[kk].x[0];
                        ytemp = atoms[kk].x[1];
                        ztemp = atoms[kk].x[2];
                        for (l = 0; l < 68; l++) {
                            n = kk + NeighborOffset[l];
                            if (atoms[n].type == -1) continue;
                            delx = xtemp - atoms[n].x[0];
                            dely = ytemp - atoms[n].x[1];
                            delz = ztemp - atoms[n].x[2];
                            dist2 = delx * delx + dely * dely + delz * delz;
                            if (dist2 < (cutoffRadius * cutoffRadius)) {
                                r = sqrt(dist2);
                                p = r * invDx + 1.0;
                                m = (int) p;
                                m = max(1, min(m, (phi_n - 1)));
                                p -= m;
                                p = min(p, 1.0);

                                _spline[2] = (((_spline_[m - 2] - _spline_[m + 2]) +
                                               8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) * invDx;
                                _spline[1] = 2.0 * (3.0 * (_spline_[m + 1] - _spline_[m]) - 2.0 * (((_spline_[m - 2] -
                                                                                                     _spline_[m + 2]) +
                                                                                                    8.0 *
                                                                                                    (_spline_[m + 1] -
                                                                                                     _spline_[m - 1])) /
                                                                                                   12.0) -
                                                    (((_spline_[m - 1] - _spline_[m + 3]) +
                                                      8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0)) * invDx;
                                _spline[0] = 3.0 * ((((_spline_[m - 2] - _spline_[m + 2]) +
                                                      8.0 * (_spline_[m + 1] - _spline_[m - 1])) / 12.0) +
                                                    (((_spline_[m - 1] - _spline_[m + 3]) +
                                                      8.0 * (_spline_[m + 2] - _spline_[m])) / 12.0) -
                                                    2.0 * (_spline_[m + 1] - _spline_[m])) * invDx;

                                dRho = (_spline[0] * p + _spline[1]) * p + _spline[2];
                                recip = 1.0 / r;
                                psip = (atoms[kk].df + atoms[n].df) * dRho;
                                fpair = -psip * recip;

                                atoms[kk].f[0] += delx * fpair;
                                atoms[kk].f[1] += dely * fpair;
                                atoms[kk].f[2] += delz * fpair;
                                atoms[n].f[0] -= delx * fpair;
                                atoms[n].f[1] -= dely * fpair;
                                atoms[n].f[2] -= delz * fpair;
                            }
                        }
                    }
                    kk++;
                }
                put_reply = 0;
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k);
                // athread_put(PE_MODE, &atoms[0], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch1 = dma_p2p(&atoms[0], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);                
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 1);
                // athread_put(PE_MODE, &atoms[nxy_per_step], &atomlist[kk] + 2, nxy_per_step * 104, (void *) &put_reply,
                //             (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch2 = dma_p2p(&atoms[nxy_per_step], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);                  
                kk = IndexOf3DIndex(nx_per_step_cal * (m_step % n_block), j - 2, k + 2);
                // athread_put(PE_MODE, &atoms[nxy_per_step * 2], &atomlist[kk] + 2, nxy_per_step * 104,
                //             (void *) &put_reply, (nghostx - nx_per_step) * 104, nx_per_step * 104);
                ch3 = dma_p2p(&atoms[nxy_per_step * 2], 1, nxy_per_step * 104, 0,
                              &atomlist[kk] + 2, ny_per_step, nx_per_step * 104, (nghostx - nx_per_step) * 104,
                             false, 0);         
                dma_wait(ch1);
                dma_wait(ch2);
                dma_wait(ch3);        
                // while (put_reply != 3);
            }
        }
    }
}

int BeginOfArea(int natom, int start, int subindex) {
    return start + subindex * block_size;
}

int EndOfArea(int natom, int start, int subindex) {
    if (subindex < Nvalue - bigblocknum) {
        return start + (subindex + 1) * block_size;
    }
    return start + natom;
}

long int IndexOf3DIndex(long int xIndex, long int yIndex, long int zIndex) {
    return (zIndex * nghosty + yIndex) * nghostx + xIndex;
}

void rtc(unsigned long *counter) {
    unsigned long rpcc;
    asm volatile("rcsr %0,4":"=r"(rpcc));
    *counter = rpcc;
}
